\documentclass[a4paper]{article}
\usepackage[affil-it]{authblk}
\usepackage{amssymb}
\usepackage{verbatim}
\usepackage{indentfirst}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{geometry}
\usepackage{xeCJK}
\usepackage{float}
\usepackage{tikz}
\usepackage{listings}
\usepackage{graphicx} % Required for figures
\usepackage{subcaption} % Required for subfigures
\usepackage{tocloft} % Required for table of contents
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

\begin{document}
% =================================================

\title{Numerical Analysis Programming Project}

\author{Zhou Chuandi 3220101409}
\affil{Qiangji Mathematics 2201}

\date{\today}

\maketitle

\begin{abstract}
This project explores the spline interpolation.

The corresponding codes of the problems in the textbook are named as \texttt{A.cpp}, \texttt{CD.cpp}, \texttt{E.cpp}, and \texttt{F.cpp}.

For the project assignment requirements, I use T1 ... T7 to refer to them.

In the bonus part, I \textbf{designed more function examples for testing}, completed together with T4. Additionally, I implemented an \textbf{algorithm to determine whether a curve is self-intersecting} (included in \texttt{CurveFitting.hpp}) and tested the curves in problem \(E\) and two other self-intersecting curves in \texttt{SP.cpp}.

The order of the sections in the report may differ from the order in the project requirements provided by the teaching assistant. I have arranged them according to my own design logic.

Happy New Year......

\end{abstract}

\tableofcontents

\section{Class Design and Inheritance}
I designed three header files \texttt{Function.hpp}, \texttt{Spline.hpp}, and \texttt{CurveFitting.hpp} and used Doxygen to generate design documentation. Go to \texttt{./Doxydocs/html/index.html}.

\section{Design of Bspline class and T5: Draw Bsplines of any degrees}
The \texttt{Spline} class contains the knots and the degree. And I would like its subclass Bspline to contain the coefficients to generate the Bspline.

T5 is ''drawing Bsplines for any given \(n\) by given knots \([t_1, \ldots, t_N]\) and coefficients \(a_i\)".
\[
S(t) = \sum_{i=2-n}^{N} a_i B_i^n(t) \in S_{n}^{n-1}(t_1, \ldots, t_n), \quad t \in [t_1, t_N]
\]

However, in fact, \(B_i^n(t)\) is determined by knots \([t_{i-1}, \ldots, t_{i+n}]\). So we need to add one knot before \(t_1\) and \(n\) knots after \(t_N\). The different choice of these extra knots will lead to different results. (It's theoretically true and I've tested as well.) So I will give knots \(t_0 \cdots t_{N+n}\) and coefficients \(a_1 \cdots a_N\) as input.

In \texttt{T5.cpp} I drew degree \(2,4,5,6,7\) Bsplines with given coefficients and knots. The results are as shown in \ref{fig:T51}.

\begin{figure}[!htbp]
  \centering
  % 第一张图片
  \begin{subfigure}{0.16\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/T5_Degree2.png}
    \caption{n=2}
    \label{fig:T511}
  \end{subfigure}
  \begin{subfigure}{0.16\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/T5_Degree4.png}
    \caption{n=4}
    \label{fig:T512}
  \end{subfigure}
  \begin{subfigure}{0.16\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/T5_Degree5.png}
    \caption{n=5}
    \label{fig:T513}
  \end{subfigure}
  \begin{subfigure}{0.16\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/T5_Degree6.png}
    \caption{n=6}
    \label{fig:T514}
  \end{subfigure}
  \begin{subfigure}{0.16\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/T5_Degree7.png}
    \caption{n=7}
    \label{fig:T515}
  \end{subfigure}
  \caption{Draw Bsplines of any degrees by given knots and coefficients}
  \label{fig:T51}
\end{figure}

\section{T1: Linear Splines}

Let \(x_i\) be the function value at \(f_i\).

\subsection{PPForm Linear Splines}

LinearPPFormSpline directly connects \((x_i,f_i)\) with line segments.

\subsection{Linear BSplines}

LinearBSpline directly sets the coefficients of the B base function to \(f_i\).

They will be tested for comparison later.

\section{T2: PPForm Cubic Splines}

For \(n\) knots \(x_0\) to \(x_{n-1}\) (\(i = 0\) to \( n-1 \) is for better coding), we take \(M_i=S''(x_i)\) as the unknowns, and we have equations to solve them.

The formula for the spline on \([x_i,x_{i+1}]\) is (using integration):
\[ S(x) = M_i \frac{(x_{i+1} - x)^3}{6h_i} + M_{i+1} \frac{(x - x_i)^3}{6h_i} + \left(y_i - \frac{M_i h_i^2}{6}\right) \frac{x_{i+1} - x}{h_i} + \left(y_{i+1} - \frac{M_{i+1} h_i^2}{6}\right) \frac{x - x_i}{h_i} \]

Where \[h_i = x_{i+1} - x_i\]

For \( i = 1, \ldots, n-2 \), we use some formulas to obtain that:
\[\lambda_i M_{i-1}+2M_i+\mu_i M_{i+1}=d_i\]

Where \[ d_i = 6f[x_{i-1},x_i,x_{i+1}]=\frac{6}{h_i + h_{i-1}} \left( \frac{y_{i+1} - y_{i}}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}} \right) \]

\[ \lambda_i = \frac{h_{i-1}}{h_i + h_{i-1}} ,\quad \mu_i = \frac{h_i}{h_i + h_{i-1}} \]

\(\mu_0\), \(\lambda_{n-1} \),\(d_0\), \(d_{n-1} \) are determined by different boundary conditions.

\subsection{Complete Cubic PPform}

We have \[ S'(x_0) = f'_0, \quad S'(x_{n-1}) = f'_{n-1} \]

Using some formulas, we obtain that \[ 2M_0 + M_1 = 6[x_0,x_0,x_1],\quad M_{n-2} + 2M_{n-1} = 6[x_{n-2},x_{n-1},x_{n-1}] \]

And we are going to solve:
\[
\begin{pmatrix}
2 & \mu_0 & & & & \\
\lambda_1 & 2 & \mu_1 & & & \\
 & \ddots & \ddots & \ddots & & \\
 & & \lambda_{n-3} & 2 & \mu_{n-3} & \\
 & & & \lambda_{n-2} & 2 & \mu_{n-2} \\
 & & & & \lambda_{n-1} & 2
\end{pmatrix}
\begin{pmatrix}
M_0 \\
M_1 \\
\vdots \\
M_{n-3} \\
M_{n-2} \\
M_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
d_0 \\
d_1 \\
\vdots \\
d_{n-3} \\
d_{n-2} \\
d_{n-1}
\end{pmatrix}
\]

Where \[\mu_0=1, d_0=\frac{6}{h_0}\left(\frac{y_1-y_0}{h_0}-f'_0\right) \]
\[ \lambda_{n-1}=1, d_{n-1}=\frac{6}{h_{n-2}}\left(f'_{n-1}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\right) \]

\subsection{Natural Cubic PPform (or Specified second derivative boundaries)}

We have \[ S''(x_0) = L, \quad S''(x_{n-1}) = R \]

And we are going to solve:
\[
\begin{pmatrix}
1 & & & & & \\
\lambda_1 & 2 & \mu_1 & & & \\
 & \ddots & \ddots & \ddots & & \\
 & & \lambda_{n-3} & 2 & \mu_{n-3} & \\
 & & & \lambda_{n-2} & 2 & \mu_{n-2} \\
 & & & & & 1
\end{pmatrix}
\begin{pmatrix}
M_0 \\
M_1 \\
\vdots \\
M_{n-3} \\
M_{n-2} \\
M_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
L \\
d_1 \\
\vdots \\
d_{n-3} \\
d_{n-2} \\
R
\end{pmatrix}
\]

\subsection{Periodic Cubic PPform}
A periodic cubic spline is used when the function is periodic, and the spline must also be periodic.

We have \(M_0 = M_{n-1}\).

Using some formulas, we obtain that \[\mu_{n-1}M_1+\lambda_{n-1}M_{n-2} + 2M_{n-1} = d_{n-1}\]

Where

\[\mu_{n-1}=\frac{h_0}{h_{n-2}+h_0},\quad \lambda_{n-1} =1 - \mu_{n-1},\quad d_{n-1}=\frac{6}{h_0+h_{n-2}}\left(\frac{y_1-y_0}{h_0}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\right)\]

And we are going to solve:
\[
\begin{pmatrix}
 1 & & & & & & -1 \\
 \lambda_1 & 2 & \mu_1 & & & & \\
 & \mu_2 & 2 & \mu_2 & & & \\
 & & \ddots & \ddots & \ddots & & \\
 & & & \lambda_{n-3} & 2 & \mu_{n-3} & \\
 & & & & \lambda_{n-2} & 2 & \mu_{n-2} \\
 & \mu_{n-1} & & & & \lambda_{n-1} & 2
\end{pmatrix}
\begin{pmatrix}
M_0 \\
M_1 \\
M_2 \\
\vdots \\
M_{n-3} \\
M_{n-2} \\
M_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
0 \\
d_1 \\
d_2 \\
\vdots \\
d_{n-3} \\
d_{n-2} \\
d_{n-1}
\end{pmatrix}
\]

\section{T3: Cubic Bsplines}

By \textbf{Theorem 3.49} Bspline \(s(x)\) of degree \(n\) over knots \(x_1 \cdots x_N\) can be expressed as:

\[s(x) = \sum_{i=2-n}^N a_i B_i^n(x) \]

\textbf{Remark: If we have \(N\) knots associated with function values, we have \(N+n-1\) knots associated with coefficients, and \(N+2n\) knots in total.}

So here we need to add \(n\) virtual knots before \(x_1\) and \(n\) virtual knots after \(x_N\). Only \(n-1\) virtual knots before \(x_1\) have non-zero coefficients.

Here \(n=3\), and since we are using \texttt{C++}, subscripts start from 0. So I use \[s(x) = \sum_{i=0}^{N+1} a_i B_{i+1}^3(x) \]

Where \(x_3 \cdots x_{N+2}\) are the original \(N\) knots given and \(f_3 \cdots f_{N+2}\) are the function values at these knots.

Remember to get the unique cubic Bspline, we need to add \(1\) knot before and \(3\) knots after. The coefficients change while the choice of extra knots changes and so the final results are the same. I just set the gap between them to be \(x_4-x_3\).

B-spline \(B_i^3(x)\) is not zero on knots \(x_{i},\ x_{i+1},\ x_{i+2}\) and zero on all other knots. So,
\[ a_{i-3}B_{i-2}^3(x_i) + a_{i-2}B_{i-1}^3(x_i) + a_{i-1}B_i^3(x_i) = f_i, \quad i = 3, 4, \dots, N+2 \]

Now we have \(N\) equations and we need two more.

\subsection{Complete Cubic BSpline}

We have \[ s'(x_3) = f'(x_3), s'(x_{N+2})=f'(x_{N+2})\]

Take derivative Bspline and use some formulas, we obtain that:

\[-3c_1 a_0 + 3 (c_1-c_2) a_1 + 3 c_2 a_2 = f'(x_3)\]
\[-3 c_{3} a_{N-1} + 3 (c_{3}-c_{4}) a_N + 3 c_{4} a_{N+1} = f'(x_{N+2})\]

Where \[c_1 = \frac{x_4-x_3}{(x_4-x_2)(x_4-x_1)},\quad c_2 = \frac{x_3-x_2}{(x_5-x_2)(x_4-x_2)} \]
\[ c_3 = \frac{x_{N+3}-x_{N+2}}{(x_{N+3}-x_{N+1})(x_{N+3}-x_N)},\quad c_4 = \frac{x_{N+2}-x_{N+1}}{(x_{N+4}-x_{N+2})(x_{N+4}-x_{N+1})}\]

And we are going to solve:
\[
\begin{pmatrix}
-3c_1 & 3c_1-3c_2 & 3c_2 & & & & \\
 & & & & -3c_3 & 3c_3-3c_4 & 3c_4 \\
 B_1^3(x_3)&  B_2^3(x_3) &  B_3^3(x_3) & & & &\\
 & B_2^3(x_4)& B_3^3(x_4) & B_4^3(x_4) & & &\\
 & & \vdots & \vdots & \vdots & & \\
 & & & B_{N-1}^3(x_{N+1}) & B_{N}^3(x_{N+1}) & B_{N+1}^3(x_{N+1}) & \\
 & & & & B_{N}^3(x_{N+2}) & B_{N+1}^3(x_{N+2}) & B_{N+2}^3(x_{N+2})
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
a_3 \\
\vdots \\
a_{N} \\
a_{N+1}
\end{pmatrix}
=
\begin{pmatrix}
f'(x_3)  \\
f'(x_{N+2})  \\
f_3 \\
f_4 \\
\vdots \\
f_{N+1} \\
f_{N+2}
\end{pmatrix}
\]

\subsection{Natural Cubic BSpline}

We have \[s''(x_3) = 0, s''(x_{N+2}) = 0 \]

Take derivative Bspline and use some formulas, we obtain that:

\[ 6 d_2 a_2 - 6 (d_1+d_2) a_1 + 6 d_1 a_0 = 0 \]
\[ 6 d_3 a_N - 6 (d_3+d_4) a_{N+1} + 6 d_4 a_{N+2} = 0\]

Where

\[d_1 = \frac{1}{(x_4-x_2)(x_4-x_1)},\quad d_2 = \frac{1}{(x_5-x_2)(x_4-x_2)} \]
\[ d_3 = \frac{1}{(x_{N+3}-x_{N+1})(x_{N+3}-x_N)},\quad d_4 = \frac{1}{(x_{N+4}-x_{N+2})(x_{N+4}-x_{N+1})} \]

And we are going to solve:
\[
\begin{pmatrix}
 6d_1 & -6d_1-6d_2 & 6d_2 & & & & \\
 & & & & 6d_3 & -6d_3-6d_4 & 6d_4 \\
 B_1^3(x_3)&  B_2^3(x_3) &  B_3^3(x_3) & & & &\\
 & B_2^3(x_4)& B_3^3(x_4) & B_4^3(x_4) & & &\\
 & & \vdots & \vdots & \vdots & & \\
 & & & B_{N-1}^3(x_{N+1}) & B_{N}^3(x_{N+1}) & B_{N+1}^3(x_{N+1}) & \\
 & & & & B_{N}^3(x_{N+2}) & B_{N+1}^3(x_{N+2}) & B_{N+2}^3(x_{N+2})
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
a_3 \\
\vdots \\
a_{N} \\
a_{N+1}
\end{pmatrix}
=
\begin{pmatrix}
0  \\
0  \\
f_3 \\
f_4 \\
\vdots \\
f_{N+1} \\
f_{N+2}
\end{pmatrix}
\]

\subsection{Periodic Cubic BSpline}

We have \[s'(x_1) = s'(x_N), s''(x_1) = s''(x_N)\]

Take derivative Bspline and use some formulas, we obtain that:
\[
-3c_1 a_0 + 3 (c_1-c_2) a_1 + 3 c_2 a_2 = -3 c_3 a_{N-1} + 3 (c_{3}-c_{4}) a_N + 3 c_{4} a_{N+1} \]
\[ 6 d_1 a_0 - 6 (d_1+d_2) a_1 + 6 d_2 a_2 = 6 d_3 a_{N+2} - 6 (d_3+d_4) a_{N+1} + 6 d_4 a_N \]

\[
\begin{pmatrix}
 6d_1 & -6d_1-6d_2 & 6d_2 & & -6d_3 & 6d_3+6d_4 & -6d_4 \\
 -3c_1 & 3c_1-3c_2 & 3c_2 & & 3c_3 & -3c_3+3c_4 & -3c_4 \\
 B_1^3(x_3)&  B_2^3(x_3) &  B_3^3(x_3) & & & &\\
 & B_2^3(x_4)& B_3^3(x_4) & B_4^3(x_4) & & &\\
 & & \vdots & \vdots & \vdots & & \\
 & & & B_{N-1}^3(x_{N+1}) & B_{N}^3(x_{N+1}) & B_{N+1}^3(x_{N+1}) & \\
 & & & & B_{N}^3(x_{N+2}) & B_{N+1}^3(x_{N+2}) & B_{N+2}^3(x_{N+2})
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
a_3 \\
\vdots \\
a_{N} \\
a_{N+1}
\end{pmatrix}
=
\begin{pmatrix}
0  \\
0  \\
f_3 \\
f_4 \\
\vdots \\
f_{N+1} \\
f_{N+2}
\end{pmatrix}
\]

\section{Problems A, C, D, F in the textbook}

\subsection{Problem A}

In \texttt{A.cpp}, I used linear PPForm(for comparison and testing) and complete cubic PPForm to interpolate the given data points and drew the graph.

(To save on images, this problem only uses complete cubic PPForm. But \texttt{T4.cpp} contains the same function as this problem. It has confirmed that the fitting effects with either PPForm or Bspline with three boundary conditions are very close.)

For the interval \([-1,1]\) with evenly spaced \(N\) nodes where \( n=  6, 11, 21, 41, 81\). The results are as \ref{fig:A1}.

\begin{figure}[!htbp]
  \centering
  % 第一张图片
  \begin{subfigure}{0.4\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/A_Cubic.png}
    \caption{Cubic PPForm}
    \label{fig:A1_1}
  \end{subfigure}
  % 第二张图片
  \begin{subfigure}{0.4\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/A_Linear.png}
    \caption{Linear PPForm}
    \label{fig:A1_2}
  \end{subfigure}
  \caption{Interpolation of Problem A against the exact function}
  \label{fig:A1}
\end{figure}

From \ref{fig:A1}, we can see that the cubic PPForm has a better fitting effect than the linear PPForm and spline interpolation is \textbf{free from the Runge phenomenon}.

The convergence rate \(k\) Satisfies that \(errors = O(N^{-k})\), so \(k\) is approximate to the tangent of the curve in the log-log plot of errors and \(N\).

From \ref{fig:A3}, we can approximate that the convergence rate of cubic splines is \(4\), and the convergence rate of linear PPForm is \(2\).

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.5\textwidth]{images/A_Convergences.png}
  \caption{Convergence Rates of Problem A}
  \label{fig:A3}
\end{figure}

\subsection{Problem C}

We were not required to do quadratic splines, so here I compared the \textbf{cubic} and the \textbf{linear} \textbf{Bspline} in \texttt{CD.cpp}.

For \( f(x)=\frac{1}{1+x^2} \) on \([-5,5]\) with \( t_i = -5, -4, \cdots, 5\), the results are as shown in \ref{fig:C1}.

\begin{figure}[!htbp]
  \centering
  \begin{subfigure}{0.4\textwidth}
    \includegraphics[width=\textwidth]{images/C.png}
    \caption{Interpolation of C against the exact function}
    \label{fig:C1}
  \end{subfigure}
  \begin{subfigure}{0.4\textwidth}
    \includegraphics[width=\textwidth]{images/D.png}
    \caption{Errors for Problem D}
    \label{fig:D1}
  \end{subfigure}
  \caption{Comparison of cubic and linear Bspline interpolations for Problems C and D}
  \label{fig:CD}
\end{figure}

\subsection{Problem D}

Problem D uses the same function as problem C so I put them both in \texttt{CD.cpp}. The results are as shown in \ref{fig:D1}.

Form \texttt{output/D.txt}, we get the errors in Table \ref{tab:D2}. 

\begin{table}[!htbp]
  \centering
  \begin{tabular}{|c|c|c|}
  \hline
  \textbf{X} & \textbf{Cubic Error} & \textbf{Linear Error} \\ \hline
  -3.5       & 0.00394007          & 0.000679392         \\ \hline
  -3         & 0                   & 1.38778e-17         \\ \hline
  -0.5       & 0.05                & 0.0205291           \\ \hline
  0          & 0                   & 0                   \\ \hline
  0.5        & 0.05                & 0.0205288           \\ \hline
  3          & 0                   & 0                   \\ \hline
  3.5        & 0.00394007          & 0.000669567         \\ \hline
  \end{tabular}
  \caption{Error Comparison between Cubic and Linear Interpolations}
  \label{tab:D2}
\end{table}

We can see the cubic Bspline has a better fitting effect than the linear Bspline. 

The errors at \(-3,0,3\) are close to machine precision. This is beacause they are the knots that generate the Bspline and the function values must be the same.

\subsection{Problem F}

The picture for \(n=1\) is shown in \ref{fig:F1} and the picture for \(n=2\) is shown in \ref{fig:F2}.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.7\textwidth]{images/F1.png}
  \caption{Divided Differences of n=1 in Problem F}
  \label{fig:F1}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.9\textwidth]{images/F2.png}
  \caption{Divided Differences of n=2 in Problem F}
  \label{fig:F2}
\end{figure}

\section{Bonus: More Function Examples}

I tested the following functions in \texttt{T4.cpp}: \(F1(x) = \frac{1}{1 + 25x^2}\), \(F2(x) = \sin(2\pi x)\), \(F3(x) = e^{-x^2}\) with \textbf{unevenly spaced knots}.
From the results in \ref{fig:T4_Collection}, we can see that:
\begin{enumerate}
  \item The fitting effects of both formats of cubic spline is much better than the linear spline.
  \item The fitting effects of the cubic spline with different boundary conditions are very close.
\end{enumerate}

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F1_Linear.png}
        \caption{$F_1(x) = \frac{1}{1 + 25x^2}, \text{Linear}$}
        \label{fig:T4F1L}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F1_Natural.png}
        \caption{$F_1(x) = \frac{1}{1 + 25x^2}, \text{Natural}$}
        \label{fig:T4F1N}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F1_Periodic.png}
        \caption{$F_1(x) = \frac{1}{1 + 25x^2}, \text{Periodic}$}
        \label{fig:T4F1P}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F1_Complete.png}
        \caption{$F_1(x) = \frac{1}{1 + 25x^2}, \text{Complete}$}
        \label{fig:T4F1C}
    \end{subfigure}
    
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F2_Linear.png}
        \caption{$F_2(x) = \sin(2\pi x), \text{Linear}$}
        \label{fig:T4F2L}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F2_Natural.png}
        \caption{$F_2(x) = \sin(2\pi x), \text{Natural}$}
        \label{fig:T4F2N}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F2_Periodic.png}
        \caption{$F_2(x) = \sin(2\pi x), \text{Periodic}$}
        \label{fig:T4F2P}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F2_Complete.png}
        \caption{$F_2(x) = \sin(2\pi x), \text{Comp-}$}
        \label{fig:T4F2C}
    \end{subfigure}
    
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F3_Linear.png}
        \caption{$F_3(x) = e^{-x^2}, \text{Linear}$}
        \label{fig:T4F3L}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F3_Natural.png}
        \caption{$F_3(x) = e^{-x^2}, \text{Natural}$}
        \label{fig:T4F3N}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/T4_F3_Periodic.png}
        \caption{$F_3(x) = e^{-x^2}, \text{Periodic}$}
        \label{fig:T4F3P}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
      \includegraphics[width=\linewidth]{images/T4_F3_Periodic.png}
      \caption{$F_3(x) = e^{-x^2}, \text{Complete}$}
      \label{fig:T4F3C}
    \end{subfigure}
    \caption{Collection of T4 Curves with Different Spline Types}
    \label{fig:T4_Collection}
\end{figure}

\section{T4: Equivalence Check of Cubic Splines}

From the results in \ref{fig:T4_Collection}, we can see that when using the same interpolation points and the same boundary conditions, the curves obtained from both the ppForm and BSpline formats are very close.

I also tested the equivalence of the two formats in \texttt{T4.cpp} by the differences between the values of two formats. The results are in \texttt{output/TestEquivalence.txt} and the maximum difference is close to machine precision, which implies that \textbf{The results generated by Cubic PPForm and Cubic Bspline are the same when using the same proper interpolation points and the same boundary conditions.}

\section{T6: Curve Fitting, and Problem E}

We have heart-shaped curves \(R1\) and Spiral \(R2\) with their parametric equations. We fit the curves with Bezier curves and cubic splines. 

In this problem, \(x(t)\) and \(y(t)\) are separately interpolated by cubic splines. Since the first derivatives of \(y(t)\) for \(R1\) do not exist at the endpoints, I chose to use a cubic natural spline (using PPForm).

Aldo I compared the different number \((N=10,N=40,N=160)\) of cumulative chordal length and equal length knots for both Bezier and cubic-spline curve fitting.

Whether using chordal length is a parameter in curvefitting classes.

From the figures in \ref{fig:EComparison}, we can see that :
\begin{enumerate}
  \item Larger \(N\) leads to a better fitting effect. For a large \(N\), the fitting effects of both the cubic spline and the Bezier curve, when using either type of nodes, are equally good.
  \item Cumulative chordal length nodes may provide a more adaptive representation of the curve, allocating more nodes in regions where the curve changes more rapidly.
  \item From figures with \(N=10\) and \(N=40\), We can see that the cubic spline has a better fitting effect when using cumulative chordal length nodes.
  \item Cubic spline interpolation ensures better smoothness. And it is more stable, as shown in \ref{fig:ER2Chordal10}
\end{enumerate}

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Chordal_10.png}
        \caption{R1, Chordal length, N=10}
        \label{fig:ER1Chordal10}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Equal_10.png}
        \caption{R1, Equal length, N=10}
        \label{fig:ER1Equal10}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Chordal_10.png}
        \caption{R2, Chordal length, N=10}
        \label{fig:ER2Chordal10}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Equal_10.png}
        \caption{R2, Equal length, N=10}
        \label{fig:ER2Equal10}
    \end{subfigure}

    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Chordal_40.png}
        \caption{R1, Chordal length, N=40}
        \label{fig:ER1Chordal40}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Equal_40.png}
        \caption{R1, Equal length, N=40}
        \label{fig:ER1Equal40}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Chordal_40.png}
        \caption{R2, Chordal length, N=40}
        \label{fig:ER2Chordal40}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Equal_40.png}
        \caption{R2, Equal length, N=40}
        \label{fig:ER2Equal40}
    \end{subfigure}

    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Chordal_160.png}
        \caption{R1, Chordal length, N=160}
        \label{fig:ER1Chordal160}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R1_Equal_160.png}
        \caption{R1, Equal length, N=160}
        \label{fig:ER1Equal160}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Chordal_160.png}
        \caption{R2, Chordal length, N=160}
        \label{fig:ER2Chordal160}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\linewidth]{images/E_R2_Equal_160.png}
        \caption{R2, Equal length, N=160}
        \label{fig:ER2Equal160}
    \end{subfigure}

    \caption{Comparison of Bezier and Cubic Spline with different nodes and methods}
    \label{fig:EComparison}
\end{figure}

\section{Bonus: Curve Self-Intersection Check Algorithm}

The algorithm to determine whether a curve is self-intersecting is implemented in \texttt{CurveFitting.hpp}. The basic idea is to discretize the curve into a number of segments and then check if any two segments intersect.

\begin{enumerate}
  \item Divide the curve into \(n\) segments by chordal length.
  \item For each pair of segments that are not adjacent, check if they are nearby to reduce computational complexity. (I used L1 distance and a threshold of \(3(|\Delta x| + |\Delta y|)\)).
  \item If they are nearby, check if two segments intersect by basic computational geometry. If they do, return true. 
  \item If no two segments intersect, return false.
  \item The time complexity is \(O(n^2)\).
\end{enumerate}

The parametric equations for the four curves are given by:
\begin{align*}
\begin{array}{l}
R1: \begin{cases}
x(t) = \sqrt{3} \cos(t) \\
y(t) = \frac{2}{3} \sqrt{|\sqrt{3} \cos(t)|} + \frac{2}{3} \sqrt{3} \sin(t)
\end{cases} \quad 
R2: \begin{cases}
x(t) = \sin(t) + t \cos(t) \\
y(t) = \cos(t) - t \sin(t)
\end{cases} \\
R3: \begin{cases}
x(t) = \cos(3t) \\
y(t) = \sin(4t)
\end{cases} \quad 
R4: \begin{cases}
x(t) = \cos(t) \\
y(t) = \sin(2t)
\end{cases}
\end{array}
\end{align*}

R1 and R2 are examples used in curve fitting and we know they are not self-intersecting. \(R3\) and \(R4\) as shown in \ref{fig:Geo3} and \ref{fig:Geo4} (From \href{https://www.geogebra.org/calculator}{GeoGebra}) are self-intersecting.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.7\textwidth]{Geo3.png}
  \caption{\(R3\)}
  \label{fig:Geo3}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.7\textwidth]{Geo4.png}
  \caption{\(R4\)}
  \label{fig:Geo4}
\end{figure}

The results of the algorithm are as follows. Correct.

\begin{lstlisting}
R1 curve is not self-intersecting

R2 curve is not self-intersecting

R3 curve is self-intersecting

R4 curve is self-intersecting
\end{lstlisting}



\section*{Acknowledgements}
1. Teacher Wang Heyu, TA Cao Shaozhen and Fan Rui, Suspicious person Hu Shuang.

2. Numerical Analysis textbooks and online resources such as \href{https://math.ecnu.edu.cn/~sfzhu/course/NumerAnal/Interpolation5.pdf}{Lecture notes from ECNU}.

3. My roommates for discussing our ideas and bugs.

4. \href{https://kimi.moonshot.cn}{Kimi AI} for additional coding assistance.

5. \href{https://chat.openai.com}{ChatGPT} (GPT-4) for code completion and debugging.

\end{document}